MODULE sefmm_interface

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: fefmm_get_J_matrix

   TYPE(fmm_basis), SAVE :: basis

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_get_J_matrix(nBas,nocc,dens,cMO,fockAO)

      USE fmm_scheme_builder, ONLY: fmm_get_scheme
      USE sefmm_cholesky_mo,  ONLY: sefmm_get_cholesky_MOs
      USE sefmm_density,      ONLY: sefmm_init_density,         &
                                    sefmm_free_density
      USE sefmm_driver,       ONLY: sefmm_get_FF_J_matrix,      &
                                    sefmm_get_NF_J_matrix,      &
                                    sefmm_get_GTO_FF_J_matrix

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)    :: nBas, nocc
      REAL(REALK),   INTENT(IN)    :: dens(nBas*(nBas+1)/2)
      REAL(REALK),   INTENT(IN)    :: cMO(nBas,nocc)
      REAL(REALK),   INTENT(INOUT) :: fockAO(nBas*(nBas+1)/2)

      TYPE(scheme_paras), POINTER :: scheme
      REAL(REALK)   :: chol_MO(nBas,nBas)
      REAL(REALK)   :: J_matrix(nBas*(nBas+1)/2)
      REAL(REALK)   :: FF_J_matrix(nBas,nBas)
      REAL(REALK)   :: NF_J_matrix(nBas,nBas)
      REAL(REALK)   :: sq_dens(nBas,nBas)
      REAL(REALK)   :: FF_energy, NF_energy, T0, T1
      INTEGER(INTK) :: i,j,ij,rank

!fixme GET BASIS SET INFO

      sq_dens = two*MATMUL(cMO,TRANSPOSE(cMO))
      CALL fmm_matrix_norm('sq_dens from cMO:     ',sq_dens,nBas*nBas)
!      DO j = 1, nocc
!      DO i = j, nbas 
!         print *, i,j, cMO(i,j)
!      END DO
!      END DO

      CALL sefmm_get_cholesky_MOs(nBas,dens,chol_MO,rank)
      sq_dens = two*MATMUL(chol_MO(:,1:rank),TRANSPOSE(chol_MO(:,1:rank)))
      CALL fmm_matrix_norm('sq_dens from chol_MO: ',sq_dens,nBas*nBas)
!      DO j = 1, rank
!      DO i = j, nbas 
!         print *, i,j, chol_MO(i,j)
!      END DO
!      END DO

!      ij = 0
!      DO j = 1, nBas
!         DO i = 1, j
!            ij = ij+1
!            sq_dens(i,j) = two*dens(ij)
!            sq_dens(j,i) = sq_dens(i,j)
!         END DO
!         sq_dens(j,j) = half*sq_dens(j,j)
!      END DO


      ! Build density for every box
      CALL fmm_get_scheme(scheme)
      !CALL sefmm_init_density(scheme,basis,sq_dens)
      !CALL sefmm_init_density(scheme,basis,cMO,nocc)
      CALL sefmm_init_density(scheme,basis,chol_MO(:,1:rank),rank)

      CALL CPU_TIME(T0)

      ! Get far-field J_matrix using FMM
      FF_J_matrix(:,:) = zero
      !CALL sefmm_get_GTO_FF_J_matrix(scheme,basis,sq_dens,FF_J_matrix)
      !CALL fmm_matrix_norm('FF_Jmat:',FF_J_matrix,SIZE(FF_J_matrix))
      CALL sefmm_get_FF_J_matrix(scheme,basis,FF_J_matrix,FF_energy)

      ! Get near-field J_matrix using finite elements
      NF_J_matrix(:,:) = zero
      CALL sefmm_get_NF_J_matrix(scheme,basis,NF_J_matrix,NF_energy)
      !CALL fmm_matrix_norm('NF_Jmat:',NF_J_matrix,SIZE(NF_J_matrix))

      CALL sefmm_free_density

!      ij = 0
!      FF_energy = zero
!      DO j = 1, nBas
!         DO i = 1, j
!            ij = ij+1
!            J_matrix(ij) = FF_J_matrix(i,j) + NF_J_matrix(i,j)
!            FF_energy = FF_energy + FF_J_matrix(i,j)*sq_dens(i,j)*half
!         END DO
!      END DO
!      fockAO(:) = fockAO(:) + J_matrix(:)

      CALL CPU_TIME(T1)
      CALL TIMTXT('>>> TIME USED get SE energy', T1-T0, LUPRI)
      WRITE(LUPRI,*) "-------------------------------------------------------"
      WRITE(LUPRI,*) "NF Coulomb energy    = ", NF_energy
      WRITE(LUPRI,*) "FF Coulomb energy    = ", FF_energy
      WRITE(LUPRI,*) "SEFMM Coulomb energy = ", FF_energy + NF_energy
      WRITE(LUPRI,*) "-------------------------------------------------------"

   END SUBROUTINE fefmm_get_J_matrix

!-------------------------------------------------------------------------------

END MODULE sefmm_interface
